/*
 * Topol.cpp
 *
 *  Created on: May 12, 2012
 *      Author: marchi
 */

#include "Topol.h"

namespace Topol_NS {
bool myunique(MyResidue & one, MyResidue & two){
	return (one.n==two.n);
}

class my_radius{
	const int nn;
	enum {H,C,N,O,S,P,Na,Cl,Mg,Ca};
	vector<double> radius;
	vector<string> rtype;
public:
	my_radius():nn(10){
		double tmp[]={0.03,0.175,0.165,0.155,0.176,0.19,0.12,0.202,0.105,0.12};
		string tmp2[]={"H","C","N","O","S","P","NA","CL","MG","CA"};
		radius=vector<double>(tmp,tmp+nn);
		rtype=vector<string>(tmp2,tmp2+nn);
	}
	double operator()(const string & x){
		try{
			for(int n=0;n<=Ca;n++){
				if(x == rtype[n]) {return radius[n];}
			}
			throw string("Stop Here: Can\'t find radius to assign to element type "+ x);
		}
		catch(const string & s){
			cout << s<< endl;
			exit(-1);
		}
		return -1.0;
	}
};
class op_compare{
	int seq, indx, offset;
	string sindx;
public:
	op_compare(int nseq,int nindx):seq(nseq), indx(nindx), offset(0),sindx(" ") {};
	op_compare():seq(0),indx(-1),offset(0), sindx(" "){};
	MyResidue operator()(MyResidue m){
		int n=m.n+offset;
		string ns=m.name;
		if(indx==-1){
			indx=n;
			sindx=ns;
		}
		else if(indx != n){
			seq++;
			indx=n;
			sindx=ns;
		} else if(sindx != ns){
			seq++;
			offset++;
			sindx=ns;
			indx=n+1;
		}
		m.n=seq;
		m.name=sindx;
		return m;
	}
};

class op_index{
	int seq, indx;
public:
	op_index(int nseq,int nindx):seq(nseq), indx(nindx){};
	op_index():seq(1),indx(-1){};
	int operator()(int m){
		if(!m) return 0;
		int n=m;
		if(indx==-1)indx=n;
		else if(indx != n){
			seq++;
			indx=n;
		}
		return seq;
	}
};
Topol::Topol(){
}
Topol::Topol(const vector<string> & y, bool x) {
	ExtractInfo(y,x);
}
Topol & Topol::operator()(const vector<string>& y,bool x){
	ExtractInfo(y,x);
	return *this;
}

Topol::~Topol() {
	// TODO Auto-generated destructor stub
}
void Topol::AllignGroups(vector<Dvect> & y, int & ng,vector<int>& bG,int tname){
	const double CUT0=1.3;
	const double CUT1=1.5;
	for(unsigned int m=0; m<cidx[tname].size(); m++){
		int n=cidx[tname][m] ;
		if(!atomname[n].compare(0,1,"HXXXXX",0,1) || bG[n]) continue;
		ng++;
		for(unsigned int p=0; p<cidx[tname].size(); p++){
			int o=cidx[tname][p] ;
			if(atomname[o].compare(0,1,"HXXXXX",0,1) || n == o || bG[o]) continue;
			if(y[n].Dist(y[o]) < CUT0){
				bG[n]=ng;
				bG[o]=ng;
			}
		}
	}

	for(unsigned int m=0; m<cidx[tname].size(); m++){
		int n=cidx[tname][m] ;
		if(bG[n]) continue;
		ng++;
		for(unsigned int p=0; p<cidx[tname].size(); p++){
			int o=cidx[tname][p] ;
			if(bG[o]|| n == o) continue;
			if(y[n].Dist(y[o]) < CUT1 ){
				bG[n]=ng;
				bG[o]=ng;
			}
		}
		if(bG[n]) continue;
		bG[n]=ng;
	}
}
void Topol::FindGroups(vector<Dvect> & y, int Water){
	int ng=0;
	bG.clear();
	bG=vector<int>(y.size(),0);
	for(unsigned int n=0;n<cidx.size();n++)
		if(n != Water) AllignGroups(y,ng,bG,n);

 	std::transform(bG.begin(),bG.end(),bG.begin(),op_index());
 	for(unsigned int m=0;m< cidx[Water].size();m++) bG[cidx[Water][m]]=0;
 	bGMax=*max_element(bG.begin(),bG.end());
 	vG=vector<vector<int> >(bGMax+1);
	for(unsigned int n=0;n<bG.size();n++)
		vG[bG[n]].push_back(n);
}
void Topol::ExtractInfo(const vector<string> & data,bool bRD){
	atomname.clear();
	resinfo.clear();

	MyResidue mytmp;
	vector<MyResidue> atres;
	vector<Dvect> xc;
	int nat=0;
	ResidueTypes mytype;
	map<int,vector<int> > cidx0;
 	for(size_t i=0;i<data.size();i++){
		if(data[i].find("ATOM") != string::npos || data[i].find("HETATM") != string::npos ) {
			string sub1=data[i].substr(11,5);
			string sub2=data[i].substr(17,4);
			string sub4=data[i].substr(76,2);
			std::stringstream(data[i].substr(22,4))>>mytmp.n;
			Dvect tmp;
			std::stringstream ss(data[i].substr(30,24));
			ss>> tmp;
			xc.push_back(tmp);

			sub1.erase(remove_if(sub1.begin(),sub1.end(),::isspace),sub1.end());
			sub2.erase(remove_if(sub2.begin(),sub2.end(),::isspace),sub2.end());
			sub4.erase(remove_if(sub4.begin(),sub4.end(),::isspace),sub4.end());
			int ntype=mytype.Sfind(sub2);
			cidx0[ntype].push_back(nat);
			if(sub4.empty()) {
				if(::isdigit(sub1[0])) {
					sub1=sub1.substr(1,sub1.size()-1)+sub1.substr(0,1);
				}
				sub4=sub1.substr(0,1);
				if(mytype.testIons(sub2)) sub4=sub1.substr(0,2);
			}

			mytmp.name=sub2;
			atomname.push_back(sub1);
			if(bRD) rd.push_back(my_radius()(sub4));
			else rd.push_back(1.0);
			atres.push_back(mytmp);
			nat++;
		}
	}
  	unsigned int typesize=cidx0.size();
 	cidx=vector<vector<int> >(typesize);
 	TypeNames=vector<string>(typesize);

 	map<int,vector<int> >::iterator it=cidx0.begin();
 	for(int n=0;it != cidx0.end();++it, n++){
 		cidx[n]=it->second;
 	}

 	for(int n=0;n<cidx0.size();n++)
 		TypeNames[n]=mytype.strFind(n);

 	cidx0.clear();
 	std::transform(atres.begin(),atres.end(),atres.begin(),op_compare());

 	for(unsigned int i=0;i<atres.size();i++){
 		resind.push_back(atres[i].n);
 	}

 	vector<MyResidue>::iterator ita=std::unique(atres.begin(),atres.end(),myunique);
 	atres.resize(ita-atres.begin());


 	for(unsigned int i=0;i<atres.size();i++){
 		resinfo.push_back(atres[i].name);
 	}
 	nr=atomname.size();
 	nres=atres.size();
 	cout << std::fixed << std::setw(6)<< "\n Topology extracted from pdb file. Number of atoms: "
 			<< nr << " Number of residue: " << nres << "\n\n";
 	cout << std::showpoint << std::fixed << std::left;
 	cout << std::setw(17) << "  Residue Name"<< std::setw(16) << "No. of Atoms" << std::setw(17)
 	<< "Residue Type" << endl;;
 	for(unsigned int n=0;n<TypeNames.size();n++){
 		cout << std::left << std::fixed<< "    " <<std::setw(10) <<
 				TypeNames[n] << "      " << std::setw(12) << cidx[n].size()
 				<< "   " << std::setw(15) << mytype.strFindResType(n) << endl;
 	}
 	maxtype=TypeNames.size();
 	for(;maxtype>0;){
 		if(cidx[maxtype-1].size()) break;
 		maxtype--;
 	}
 	FindGroups(xc,mytype.gWaterId());
}
} /* namespace Topol_NS */
